COMMENT ⊗   VALID 00005 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	begin "unplot" comment An algorithm to undiscretize data.
C00006 00003	Cubic spline procedure
C00011 00004	The main algorithm
C00015 00005	The driver program
C00017 ENDMK
C⊗;
begin "unplot" comment An algorithm to undiscretize data.

Given points $a↓0$, $\ldotss$, $a↓{n+1}$ (typically all integers with
$|a↓{j+1}-a↓j|≤1$, but this is not required), this algorithm finds
points $z↓0$, $\ldotss$, $z↓{n+1}$ such that $|z↓j-a↓j|≤{1\over2}$ for all
$j$, and also $z↓j=a↓j$ for $j = 0,n$. The points $z↓j$ vary smoothly, in
the sense that they lie on $m$ cubic polynomials $f↓k(j)$, where $z↓j=f↓k(j)$
when $x↓{k-1}≤j≤x↓k$. Here $0=x↓0<x↓1<\cdots<x↓m=n$, and the algorithm tends to
minimize $m$. It solves the problem in such a way that the $x↓k$ are integers,
and such that $f↓1(x↓0)=a↓0$, $f↓1(x↓1)=f↓2(x↓1)$, $\ldotss$, $f↓{m-1}(x↓{m-1})=
f↓m(x↓{m-1})$, $f↓m(x↓m)=a↓n$ at the junction points. Furthermore the
piecewise cubics vary smoothly in the sense that the derivatives satisfy
$f↑\prime↓k(x↓k)=f↑\prime↓{k+1}(x↓k)$ and $f↑{\prime\prime}↓k(x↓k)=f↑{\prime\prime}
↓{k+1}(x↓k)$ for $1≤k<m$. The first derivatives $f↑\prime↓1(0)=z↑\prime↓0$ and
$f↑\prime↓m(n)=z↑\prime↓n$ will be set to any desired values input to the algorithm.

The driver program incorporated in this version of the algorithm accepts data
as a sequence of $n+2$ characters that are either "+" or "0" or "-", specifying
differences or derivatives. For example, {\tt++00+00++00--0} means $z↑\prime↓0=+1$,
$a↓0=0$, $a↓1=1$, $a↓2=1$, $a↓3=1$, $a↓4=2$, $a↓5=2$, $a↓6=2$, $a↓7=3$, $a↓8=4$,
$a↓9=4$, $a↓{10}=4$, $a↓{11}=3$, $a↓{12}=2$, $z↑\prime↓{12}=0$, $n=12$.
;
require "⊂⊃⊂⊃" delimiters; "used for macros"
define # = ⊂;comment⊃; "used henceforth instead of quoted comments like this"
define nextline = ⊂('15&'12)⊃ # carriage-return and line-feed in print commands;
define thru = ⊂step 1 until⊃ # abbreviation for for clauses;
define DEBUGONLY = ⊂⊃ # changed to ⊂comment⊃ when not debugging;
define saf = ⊂safe⊃ # used when an array is believed to require no bounds checks;
DEBUGONLY redefine saf = ⊂⊃ # when debugging, belief turns to disbelief;
DEBUGONLY external procedure bail # the SAIL debugger in case of need;

comment The following variables are used for input and output by the algorithm;
integer n # size of input data;
define cap=100 # capacity of the arrays, must be $>n$;
real saf array a[0:cap] # given data $a↓j$;
real z0p, znp # given derivatives at the endpoints;
real saf array z[0:cap] # the answers $z↓j$;

comment The following variables are used as indices during the calculations;
integer j # runs from 0 to $n$;
integer k # runs from 0 to $m$;
comment Cubic spline procedure;

comment The following variables are used as input and output by the spline procedure;
integer m # number of pieces in piecewise cubic;
real saf array x[0:cap] # given junction points $x↓k$;
real saf array y[0:cap] # given function values $y↓k$;
real saf array c[0:cap] # auxiliary coefficients $c↓k$ in recurrence for $β↓k$;
real saf array d[0:cap] # reciprocal difference $d↓k$;
real saf array alf[0:cap] # derivative parameter $α↓k$;
real saf array bet[0:cap] # derivative parameter $β↓k$;
real saf array gam[0:cap] # auxiliary coefficients $\gamma↓k$ in recurrence for $β↓k$;
real saf array del[0:cap] # straight-line derivatives $\delta↓k$;
real saf array ddel[0:cap] # $\delta↓{k+1}-\delta↓k$;

procedure compute_splines;
begin comment This is a general procedure to compute cubic splines as defined
above, given $x[0]$, $\ldotss$, $x[m]$, $y[0]$, $\ldotss$, $y[m]$, \.{z0p} and
\.{znp}, where $y[k]$ is the desired value at $x[k]$, \.{z0p} is the desired
derivative at $x[0]$, and \.{znp} is the desired derivative at $x[m]$.
The cubic polynomial $f↓k(x)$ is written in the form $y{k-1}+(x-x↓{k-1})\delta↓k
+(x-x↓{k-1})↑2(x-x↓k}α↓kd↓k↑2+(x-x↓{k-1})(x-x↓k)↑2β↓kd↓k↑2$, where $\delta↓k=
(y↓k-y↓{k-1})/(x↓k-x↓{k-1})$ and $d↓k=1/(x↓k-x↓{k-1})$, and where $α↓k$ and $β↓k$
are parameters to be determined because of the conditions on derivatives.
We have
$$\eqalign{f↑\prime↓k(x↓{k-1})⊗=\delta↓k+β↓k,\cr f↑\prime↓k(x↓k)⊗=\delta↓k+α↓k,\cr}
\qquad\eqalign{f↑{\prime\prime}↓k(x↓{k-1})⊗=-2d↓k(2β↓k+α↓k),\cr
f↑{\prime\prime}↓k(x↓k)⊗=+2d↓k(β↓k+2α↓k).\cr}$$
Defining $\delta↓0=0$, $α↓0=\.{z0p}$, $\delta↓{m+1}=0$, and $β↓{m+1}=\.{znp}$,
we have $\delta↓k+α↓k=\delta↓{k+1}+β↓{k+1}$ for $0≤k≤m$, by matching first
derivatives. It remains to determine $β↓2$, $\ldotss$, $β↓m$ by using the
equations $d↓k(β↓k+2α↓k)+d↓{k+1}(2β↓{k+1}+α↓{k+1})=0$ for $1≤k<m$ to express
$β↓k$ in the form $c↓k-\gamma↓kβ↓{k+1}$ for increasing $k$, then to use this
recurrence for decreasing $k$.
;
del[0]←0;
for k←1 thru m do
	begin d[k]←1/(x[k]-x[k-1]); del[k]←(y[k]-y[k-1])*d[k];
	ddel[k-1]←del[k]-del[k-1];
	end;
del[m+1]←0; ddel[m]←-del[m];
c[1]←z0p-del[1]; gam[1]←0;
for k←1 thru m-1 do
	begin real denom # 1/denominator in calculation of $c↓k$ and $\gamma↓k$;
	denom←1/(2*(d[k]+d[k+1])-d[k]*gam[k]);
	gam[k+1]←d[k+1]*denom;
	c[k+1]←-(2*d[k]*ddel[k]+d[k+1]*ddel[k+1]+c[k]*d[k])*denom;
	end;
bet[m+1]←znp;
for k←m step -1 until 1 do
	begin bet[k]←c[k]-gam[k]*bet[k+1];
	alf[k]←bet[k+1]+ddel[k];
	end;
end;
comment The main algorithm;

procedure main_algorithm;
begin comment The idea of the main algorithm is to try first to solve the
problem with one cubic. If that fails, new junction points are added one at
a time at the place where the piecewise cubic now deviates most from the
input data. When a new junction point is added at $j$, the value $z↓j$ is
set to $a↓j+{1\over2}$, $a↓j$, or $a↓j-{1\over2}$, depending on whether or
not $(a↓{j-1}+a↓{j+1})/2$ is $>a↓j$, $=a↓j$, or $<a↓j$.
;
m←1; x[0]←0; x[1]←n; y[0]←0; y[1]←a[n];
while true do
	begin comment Main interation, either returns or increases $m$;
	integer maxloc # value of $j$ that maximizes $|z↓j-a↓j|$;
	real maxd # maximum of $|z↓j-a↓j|$;
	real newval # tentative value of $y↓k$ at new junction point $x↓k$;
	compute_splines;
	k←1; j←0;
	while k≤m do 
		begin comment Evaluation of the $z↓j$ with current splines;
		real jj # value at which we wish to evaluate $f↓k$;
		real alpha,beta # coefficients in $αx-β=((x-x↓{k-1})α↓k+
			(x-x↓k)β↓k)d↓k↑2$;
		alpha←(alf[k]+bet[k])*d[k]↑2;
		beta←(alf[k]*x[k-1]+bet[k]*x[k])*d[k]↑2;
		jj←j;
		z[j]←((alpha*jj-beta)*(jj-x[k])+del[k])*(jj-x[k-1])+y[k-1];
		jj←jj+1;
		z[j+1]←((alpha*jj-beta)*(jj-x[k])+del[k])*(jj-x[k-1])+y[k-1];
		jj←jj+1;
		z[j+2]←((alpha*jj-beta)*(jj-x[k])+del[k])*(jj-x[k-1])+y[k-1];
		if x[k]>j+3 then
			begin comment Tabulating a cubic;
			real f,df,ddf,dddf # differences;
			f←z[j+2] # $f(j+2)$;
			df←f-z[j+1] # $\Delta f(j+1)$;
			ddf←df-z[j+1]+z[j] # $\Delta↑2 f(j)$;
			dddf←6*alpha # $\Delta↑3 f$;
			j←j+3;
			while j<x[k] do
				begin ddf←ddf+dddf # $\Delta↑2 f(j-2)$;
				df←df+ddf # $\Delta f(j-1)$;
				f←f+df # $f(j)$;
				z[j]←f;
				j←j+1;
				end;
			end;
		j←x[k];
		k←k+1;
		end;
	maxloc←0; maxd←0;
	for j←1 thru n-1 do if abs(z[j]-a[j])>maxd then
		begin maxd←abs(z[j]-a[j]); maxloc←j;
		end;
	print(nextline,"Values with m=",m,":",nextline);
	for j←0 thru n do print(z[j]);
	if maxd≤0.5 then return;
	print(nextline,"Adding a junction point at ",maxloc,".");
	m←m+1; 	k←m;
	while x[k-1]>maxloc do
		begin x[k]←x[k-1]; y[k]←y[k-1];
		k←k-1;
		end;
	x[k]←maxloc;
	newval←(a[maxloc-1]+a[maxloc+1])/2;
	if newval<a[maxloc] then y[k]←a[maxloc]-.5
	else if newval>a[maxloc] then y[k]←a[maxloc]+.5
	else y[k]←newval;
	end;
end;
comment The driver program;

string s # given string;
integer char # current input character;

setprint("unrnd.tmp","B");
while true do
	begin print(nextline,"Input:"); s←inchwl;
	DEBUGONLY if s="b" then bail;
	if s=0 then done;
	setprint(null,"F"); print(s); setprint(null,"B") # echo the input;
	case (char←lop(s)) of begin
	["+"] z0p←+1;
	["0"] z0p←0;
	["-"] z0p←-1;
	else begin print("?"); z0p←0 end
	  end;
	n←0; a[0]←0;
	while true do
		begin case (char←lop(s)) of begin
		["+"] a[n+1]←a[n]+1;
		["0"] a[n+1]←a[n];
		["-"] a[n+1]←a[n]-1;
		else begin print("?"); a[n+1]←a[n] end
		  end;
		if s=0 then done;
		n←n+1;
		if n≥cap then
			begin print(nextline,"Table size exceeded!"); done;
			end;
		end;
	znp←a[n+1]-a[n]; z[n]←a[n];
	main_algorithm;
	end;
end